#include "structCells.hpp"
#include "structMesh.hpp"
#include "geometry.hpp"
#include "signalCodes.hpp"
#include "util.hpp"

namespace hms {

StructCells::StructCells( const StructMesh& mesh ) :
	StructCellsBase<StructCells>(mesh)
{}

auto StructCells::centroidsImp() const -> const Matrix2Xs& { return m_centroids; }
auto StructCells::areasImp    () const -> const ArrayXs&   { return m_areas; }
auto StructCells::circumfsImp () const -> const ArrayXs&   { return m_circumfs; }

auto StructCells::centroid(Index i) const -> Matrix2Xs::ConstColXpr {
	return m_centroids.col(i);
}

auto StructCells::area    (Index i) const -> scalar {
	return m_areas(i);
}

auto StructCells::circumf (Index i) const -> scalar {
	return m_circumfs[i];
}


void StructCells::storeData() {
	Index nx { mesh().nCellsX() };
	Index ny { mesh().nCellsY() };
	static constexpr Index ng {1};

	/* size is not equal to number of cells, but rather the indexing strategy */
	m_areas.resize( (nx + 2*ng ) * (ny + 2*ng ) );

	m_centroids.resize( NoChange, m_areas.size() );
	m_circumfs.resizeLike( m_areas );

	mesh().forEachDomainCell( [&]( const StructCell& cell ){
		m_areas( cell.index() ) = verticesToArea( mesh(), cell.vertices() );
		m_centroids.col( cell.index() ) = verticesToCentroid(
			mesh(), cell.vertices(), m_areas( cell.index() )
		);
		m_circumfs[ cell.index() ] =
			mesh().edges().lengths()( cell.edges() ).sum();
	} );

	mesh().forEachGhostCell( [&]( const StructCell& ghostCell ){
		/* get domain cell for its data */
		Index i { ghostCell.i() }, j { ghostCell.j() };
		// auto [iDom,jDom] { ghostCell.domainCell() };
		// StructCell domainCell { mesh().cell(iDom, jDom) };
		StructCell domainCell { ghostCell.domainCell() };
		Index iDom { domainCell.i() }, jDom { domainCell.j() };

		/* store ghost cell area */
		m_areas( ghostCell.index() ) = m_areas( domainCell.index() );

		/* get domain cell centroid */
		Matrix2Xs::ColXpr c { m_centroids.col( domainCell.index() ) };

		auto edge { ParamLine::Through(
			mesh().vertex(ghostCell.vertices()[0]),
			mesh().vertex(ghostCell.vertices()[1])
		) };

		/* mirror centroid along edge */
		m_centroids.col( ghostCell.index() ) = c + 2*(edge.projection(c)-c)*
			abs(i-iDom + j-jDom); /* either of those differences is always zero,
			* and the other one is the current ghost layer */
		/* same circumference for ghost cell as for neighbouring domain cell */
		m_circumfs[ ghostCell.index() ] = m_circumfs[ domainCell.index() ];
	} );
	
	auto setCornersToSignal = [&]( auto& vals, int nVals ){
		auto valsReshaped = vals.array().reshaped(nVals*(nx+2*ng), ny+2*ng);
		valsReshaped.bottomLeftCorner (nVals*ng,ng) = signalCode::uninitialised;
		valsReshaped.bottomRightCorner(nVals*ng,ng) = signalCode::uninitialised;
		valsReshaped.topLeftCorner    (nVals*ng,ng) = signalCode::uninitialised;
		valsReshaped.topRightCorner   (nVals*ng,ng) = signalCode::uninitialised;
	};
	setCornersToSignal( m_areas, 1 );
	setCornersToSignal( m_centroids, 2 );
	setCornersToSignal( m_circumfs, 1 );
}



}